Mini Project 1 & 2
Michaela Gustafson
10/16/2021
Spatial Data in R: Mini Project 1
The goals for this phase of your final project are:
Articulate an interesting research question based on a dataset you’d like to learn more about.
Develop a spatial database that contains potentially relevant explanatory variables that you’d like to explore in the context of that research question.
Demonstrate an understanding of the various workflow elements involved in designing and constructing a spatial database for subsequent visualization and analysis.
My Question: Where is the best place to go in New Zealand during the zombie apocalypse to avoid infection and survive? Specifically, I will look at 1. Where the most zombie attacks are happening 2. Are they happening in more urban or rural locations? 3. Highlight rural locations as areas of survival (the less people around, the less chance of encountering an infected person and becoming a zomibe) 4. Number of pharmacies in each city>>region to inform survivors of the best places to go where they will have regular access to medical supplies
I will answer this question using the following data:
- Census information (population, population density)
- Urban and rural district classifications
- Land area (sqkm)
- Locations and number of zombie attacks
library(here)## here() starts at /home/michaela/R/Mini-Project-1-michaelagustafsonlibrary(tidyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(terra)## terra version 1.4.22##
## Attaching package: 'terra'## The following object is masked from 'package:dplyr':
##
## src## The following object is masked from 'package:tidyr':
##
## extractlibrary(sf)## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(tidyverse)## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x terra::extract() masks tidyr::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x purrr::simplify() masks terra::simplify()
## x terra::src() masks dplyr::src()library(nngeo)
library(data.table)##
## Attaching package: 'data.table'## The following object is masked from 'package:purrr':
##
## transpose## The following objects are masked from 'package:terra':
##
## copy, shift## The following objects are masked from 'package:dplyr':
##
## between, first, lastlibrary(tmap)## Registered S3 methods overwritten by 'stars':
## method from
## st_bbox.SpatRaster sf
## st_crs.SpatRaster sflibrary(viridis)## Loading required package: viridisLitelibrary(rgdal)## Loading required package: sp## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## project# library(ggmap)A. Dependent Variable with Point Geometry: Zombie Attacks
Data source: kaggle.com
zom <- read.csv("zombies.csv")
head(zom); dim(zom)## zombieid zombie age sex rurality household water food medication tools
## 1 1 Human 18 Female Rural 1 0 Food Medication No tools
## 2 2 Human 18 Male Rural 3 24 Food Medication tools
## 3 3 Human 18 Male Rural 4 16 Food Medication No tools
## 4 4 Human 19 Male Rural 1 0 Food Medication tools
## 5 5 Human 19 Male Urban 1 0 Food Medication No tools
## 6 6 Human 19 Female Urban 1 0 Food Medication tools
## firstaid sanitation clothing documents
## 1 First aid supplies Sanitation Clothing <NA>
## 2 First aid supplies Sanitation Clothing <NA>
## 3 First aid supplies Sanitation Clothing <NA>
## 4 No first aid supplies Sanitation Clothing <NA>
## 5 First aid supplies Sanitation <NA> <NA>
## 6 First aid supplies Sanitation Clothing <NA>## [1] 200 14zombie.shp <- read_sf(here("New_Zealand_Zombie_Attack_Data.shp"))
plot(zombie.shp) #remember to use st_geometry to avoid plotting every column in the dataframe
B. Regions
1. New Zealand Region Data
Data source: https://datafinder.stats.govt.nz/layer/98765-regional-council-2019-clipped-generalised/
The regions of New Zealand are my ‘regions’ I will use to organize higher order variables.
# Find "Regional Council 2018 Clipped (generalised)"
# select the GeoPackage option in the "Vectors/tables" dropdown
# at https://datafinder.stats.govt.nz/data/ (requires registration)
# Save the result as:
unzip("statsnzregional-council-2019-clipped-generalised-SHP.zip")
nz.full.shp <- st_read("regional-council-2019-clipped-generalised.shp")## Reading layer `regional-council-2019-clipped-generalised' from data source
## `/home/michaela/R/Mini-Project-1-michaelagustafson/regional-council-2019-clipped-generalised.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 17 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1089970 ymin: 4747987 xmax: 2470042 ymax: 6223156
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000plot(st_geometry(nz.full.shp))
# print(object.size(nz_full), units = "Kb") # 14407.2 Kb
# nz <- st_simplify(nz_full, dTolerance = 1000)
names(nz.full.shp)## [1] "REGC2019_V" "REGC2019_1" "LAND_AREA_" "AREA_SQ_KM" "Shape_Leng"
## [6] "geometry"nz.full.shp$REGC2019_1## [1] "Northland Region" "Auckland Region"
## [3] "Waikato Region" "Bay of Plenty Region"
## [5] "Gisborne Region" "Hawke's Bay Region"
## [7] "Taranaki Region" "Manawatu-Wanganui Region"
## [9] "Wellington Region" "West Coast Region"
## [11] "Canterbury Region" "Otago Region"
## [13] "Southland Region" "Tasman Region"
## [15] "Nelson Region" "Marlborough Region"
## [17] "Area Outside Region"nz.reg.shp = filter(nz.full.shp, REGC2019_1 != "Area Outside Region") %>%
select(Region = REGC2019_1, AreaSQKM = AREA_SQ_KM, geometry)
st_is_valid(nz.reg.shp)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUEst_crs(nz.reg.shp)## Coordinate Reference System:
## User input: NZGD2000 / New Zealand Transverse Mercator 2000
## wkt:
## PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",
## BASEGEOGCRS["NZGD2000",
## DATUM["New Zealand Geodetic Datum 2000",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4167]],
## CONVERSION["New Zealand Transverse Mercator 2000",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",173,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",1600000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",10000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["New Zealand - onshore"],
## BBOX[-47.33,166.37,-34.1,178.63]],
## ID["EPSG",2193]]C. 1. Census Data
Data source: https://www.stats.govt.nz/topics/census
I will use census data to determine the population density/sqkm in areas of the regions (urban vs rural) and overall population density/sqkm of a region. This information will be best used to determine good regions for individuals to go to be in less populated areas and have a better chance of avoiding other people and potential zombie infection.
# read in census data
nz.census <- read.csv("Census_Usually_resident_population_count_and_change_by_region_2006_2013_and_2018.csv")
# filter to get most recent census population (2018)
names(nz.census)## [1] "X...Census.year" "Region" "Region.Code" "Measure"
## [5] "Value" "Value.Unit" "Value.Label" "Null.Reason"nz.census <- filter(nz.census, X...Census.year == "2018") %>%
select(Year = X...Census.year, Region, Population = Value)
nz.census$Region = paste(nz.census$Region, "Region", sep = " ")
nz.census[8, 2] = "Manawatu-Wanganui Region"You’ll want to leave yourself a few more “breadcrumbs”, I think. Anytime you are directly modifying the data (as you do above), it’s good to leave a note to remind you why you did
2. Urban vs Rural shapefile
Data source: https://datafinder.stats.govt.nz/layer/98743-urban-rural-2019-clipped-generalised/
unzip("statsnzurban-rural-2019-clipped-generalised-SHP.zip")
urbrur <- read_sf("urban-rural-2019-clipped-generalised.shp")
plot(st_geometry(urbrur))
urban <- urbrur %>%
select(Name = UR2019_V_1, Urban = IUR2019__1, AreaSQKM = AREA_SQ_KM, geometry)
unique(urban$Urban)## [1] "Inland water" "Large urban area" "Major urban area"
## [4] "Medium urban area" "Oceanic" "Rural other"
## [7] "Rural settlement" "Small urban area"# neet to change columsn to say just urban or rural
urban$Urban[grep("urban", urban$Urban)] <- "Urban"
urban$Urban[grep("Rural", urban$Urban)] <- "Rural"
urban$Urban[grep("water", urban$Urban)] <- "Other"
urban$Urban[grep("Oceanic", urban$Urban)] <- "Other"3. Pharmacy Location Data
Data source: https://fyi.org.nz/request/10183-total-number-of-licenced-pharmacies-in-nz-2019
Does this still count as tabular data even though I turned it into spatial data??
# read in pharmacy csv with location data
nz.pharm <- read.csv(here("nz.pharm.csv"))
nz.pharm <- na.omit(nz.pharm)
nz.pharm.sf <- st_as_sf(nz.pharm, coords = c("lon", "lat"), crs = 4326)
st_crs(nz.pharm.sf)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]nz.pharm.proj <- st_transform(nz.pharm.sf, crs = st_crs(nz.reg.shp))
st_crs(nz.reg.shp) == st_crs(nz.pharm.proj)## [1] TRUEplot(st_geometry(nz.reg.shp))
plot(st_geometry(nz.pharm.proj), add = TRUE)
nz.pharm.proj$type <- "Pharmacy"
#looks alrightD. Creating Database
Spatial Join
Questions to Answer: 1. What is the percentage of urban and rural areas within each region? 2. What is the number of attacks in each urban type within each region? 3. What is the population density of each region (people/sqkm)? 4. What is the number of attacks per population in each region? (attacks/people)
# make sure CRS are the same
st_crs(zombie.shp) == st_crs(urban)## [1] FALSEzombie.shp.proj <- st_transform(zombie.shp, crs = st_crs(urban))
st_crs(zombie.shp.proj) == st_crs(urban)## [1] TRUEst_crs(urban) == st_crs(nz.reg.shp)## [1] TRUE# ended up using a nearest neighbor approach because I kept getting points 800 something rows when I started with 668 so something was getting duplicated:
#zombie.urban.join <- st_join(urban, zombie.shp.proj, join = st_nn, k = 1, maxdist = 50)
#zombie.urban.region.join <- st_join(zombie.urban.join, nz.reg.shp, join = st_within)
#str(zombie.urban.region.join)
#zombie.urban.region.ng <- st_drop_geometry(zombie.urban.region.join)
#zombie.urban.region.df <- zombie.urban.region.ng %>%
# select(Urban, AreaSQKM.x, Attacks, Region)
#zombie.urban.region.df <- as.data.frame(zombie.urban.region.df)
#zombie.urban.region.df$Region <- as.factor(zombie.urban.region.df$Region)
#zombie.urban.region.df$Attacks <- replace_na(zombie.urban.region.df$Attacks, 0)
#zombie.urban.region.df$Attacks <- as.numeric(zombie.urban.region.df$Attacks)
#summary.df <- zombie.urban.region.df %>%
#group_by(Region, Urban) %>%
#summarise(., sum(Attacks), sqkm_dis_total = sum(AreaSQKM.x))
#summary.df <- rename(summary.df, Attacks = "sum(Attacks)")
#summary.df$urban_attacks <- paste(summary.df$Urban, summary.df$Attacks, sep = ";")
#summary.df.wide <- spread(summary.df, Urban, sqkm_dis_total, "0")
#summary.df.wide <- spread(summary.df.wide, urban_attacks, Attacks)
#str(summary.df.wide)
#summary.df.wide$Other <- as.numeric(summary.df.wide$Other)
#summary.df.wide$Rural <- as.numeric(summary.df.wide$Rural)
#summary.df.wide$Urban <- as.numeric(summary.df.wide$Urban)
#summary.df.wide <- summary.df.wide[!is.na(summary.df.wide$Region),]
#summary.df.wide[is.na(summary.df.wide)] <- 0
#summary.df.wide <- summary.df.wide %>%
# group_by(Region) %>%
#summarise_if(., is.numeric, sum)
#summary.df.wide <- rename(summary.df.wide, Other_sqkm = "Other")
#summary.df.wide <- rename(summary.df.wide, Rural_sqkm = "Rural")
#summary.df.wide <- rename(summary.df.wide, Urban_sqkm = "Urban")
#
#colnames(summary.df.wide)[grepl('Other;',colnames(summary.df.wide))] <- 'Other_attacks'
#colnames(summary.df.wide)[grepl('Rural;',colnames(summary.df.wide))] <- 'Rural_attacks'
#colnames(summary.df.wide)[grepl('Urban;',colnames(summary.df.wide))] <- 'Urban_attacks'
#names(summary.df.wide) <- make.unique(names(summary.df.wide))
#tidy.summary.df <- summary.df.wide %>%
#rowwise() %>%
#mutate(Rural_Attack_sum = sum(across(starts_with("Rural_attacks")), na.rm = T))
#tidy.summary.df <- tidy.summary.df %>%
#rowwise() %>%
# mutate(Other_Attack_sum = sum(across(starts_with("Other_attacks")), na.rm = T))
#tidy.summary.df <- tidy.summary.df %>%
#rowwise() %>%
#mutate(Urban_Attack_sum = sum(across(starts_with("Urban_attacks")), na.rm = T))
# taking out Other Attack sum becuase there were none (they were in an unidentified region)
#tidy.summary.df <- tidy.summary.df %>%
#select(Region, Other_sqkm, Rural_sqkm, Urban_sqkm, Rural_Attack_sum, Urban_Attack_sum)
### NEED TO FIX REGIONAL SQKM# Find the number of pharmacies per region
pharm.per.region <- st_join(nz.reg.shp, nz.pharm.proj, join = st_contains)
pharm.per.region.ng <- st_drop_geometry(pharm.per.region)
pharm.per.region.ng <- pharm.per.region.ng %>%
select(Region, Pharmacy_Nm)
# Count the number of pharmacies in each region
pharm.per.region.count <- pharm.per.region.ng %>%
group_by(Region) %>%
tally()
pharm.per.region.count <- rename(pharm.per.region.count, pharms_per_region = n)
# combine pharmacies per region with summary.df
#tidy.summary.df <- left_join(tidy.summary.df, pharm.per.region.count)Add census data to tidy summary df
# add total population to final df
#tidy.summary.df <- left_join(tidy.summary.df, nz.census)
# add total area (sqkm) to final df
#tidy.summary.df <- left_join(tidy.summary.df, nz.reg.shp)
# calculate attacks/rural area and attacks per urban area per region
#tidy.summary.df <- tidy.summary.df %>%
#mutate(attacks_per_rural = Rural_Attack_sum/Rural_sqkm,
# attacks_per_urban = Urban_Attack_sum/Urban_sqkm,
# pop_density = Population/AreaSQKM,
# perc_rural = Rural_sqkm/AreaSQKM,
# perc_urb = Urban_sqkm/AreaSQKM,
# total_attacks = Rural_Attack_sum + Urban_Attack_sum)<<<<<<< HEAD # Add elevation
rastlist1 <- list.files(path = "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/", pattern = '.tif', full.names=TRUE)
rastlist1## [1] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//aklnd_25r.tif"
## [2] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//ecape_25r.tif"
## [3] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//hksbay_25r.tif"
## [4] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//nthcape_25r.tif"
## [5] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//taranaki_25r.tif"
## [6] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//waikato_25r.tif"
## [7] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//well_25r.tif"
## [8] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//whngrei_25r.tif"rn1 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/aklnd_25r.tif")
rn2 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/ecape_25r.tif")
rn3 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/hksbay_25r.tif")
rn4 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/nthcape_25r.tif")
rn5 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/taranaki_25r.tif")
rn6 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/waikato_25r.tif")
rn7 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/well_25r.tif")
rn8 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/whngrei_25r.tif")
rastlist2 <- list.files(path = "/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/", pattern = '.tif', all.files=TRUE, full.names=FALSE)
rastlist2## [1] "chch25r.tif" "dunedin_25r.tif" "greymth_25r.tif" "invcgll_25r.tif"
## [5] "kaik_25r.tif" "mtcook_25r.tif" "nelson_25r.tif" "teanau_25r.tif"
## [9] "waitaki_25r.tif"rs1 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/chch25r.tif")
rs2 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/dunedin_25r.tif")
rs3 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/greymth_25r.tif")
rs4 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/invcgll_25r.tif")
rs5 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/kaik_25r.tif")
rs6 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/mtcook_25r.tif")
rs7 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/nelson_25r.tif")
rs8 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/teanau_25r.tif")
rs9 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/waitaki_25r.tif")
all.rast <- merge(rn1, rn2, rn3, rn4, rn5, rn6, rn7, rn8, rs1, rs2, rs3, rs4, rs5, rs6, rs7, rs8, rs9)
zombie.elev.proj <- st_transform(zombie.shp, crs = crs(all.rast))
plot(all.rast)
plot(zombie.elev.proj, add = TRUE)## Warning in plot.sf(zombie.elev.proj, add = TRUE): ignoring all but the first
## attribute
# Technically should be calculating the mean elevation or range or something from a buffer surrounding the point, but for the sake of time (because I already tried it and I'd have to mess around with changing extents and a bunch of other stuff) and to actually practice what I'm supposed to be practicing - just going to the pull the elevation from the exact point of attack:
zomb.elev.proj.vect <- vect(zombie.elev.proj)
elev.at.attack <- extract(all.rast, zomb.elev.proj.vect, na.rm = TRUE)
elev.at.attack$ID <- zombie.shp$OBJECTID
elev.at.attack <- rename(elev.at.attack, OBJECTID = ID)Trying to create new df
# add urban or rural
zombie.attacks.df <- st_join(zombie.shp.proj, urban, join = st_within)
colnames(zombie.attacks.df)## [1] "OBJECTID" "OBJECTID_1" "Day" "Attacks" "Time"
## [6] "Longitude" "Latitude" "geometry" "Name" "Urban"
## [11] "AreaSQKM"zombie.attacks.df <- zombie.attacks.df[,c(1, 4, 6:8, 10)]
# add elevation
zombie.attacks.df <- left_join(zombie.attacks.df, elev.at.attack)## Joining, by = "OBJECTID"zombie.attacks.df <- rename(zombie.attacks.df, elev = aklnd_25r)
# add region
zombie.attacks.df <- st_join(zombie.attacks.df, nz.reg.shp, join = st_within)
# add number of pharmacies per region
zombie.attacks.df <- left_join(zombie.attacks.df, pharm.per.region.count)## Joining, by = "Region"# add census data
zombie.attacks.df <- left_join(zombie.attacks.df, nz.census)## Joining, by = "Region"# mutate and calculate population density per region
zombie.attacks.df <- zombie.attacks.df %>%
mutate(pop_density = Population/AreaSQKM)
#####
nz.reg.shp.join <- left_join(nz.reg.shp, nz.census)## Joining, by = "Region"nz.reg.shp.join <- nz.reg.shp.join %>%
mutate(pop_density = Population/AreaSQKM)======= you didn’t have any raster extractions here and the analysis is not reproducible because a) you didn’t save the nz_pharm object, b) you didn’t provide all of the necessary info to recreate it. I would have liked to see a little more effort and annotation. If I knew more about what you were trying to do, I might actually be able to fix the parts of this that wont run.
save.image(here("workspaces/MiniProj_Workspace.RData"))Mini Project 2
library(tidyverse)
library(pander)##
## Attaching package: 'pander'## The following object is masked from 'package:terra':
##
## wraplibrary(sf)
library(units)## udunits database from /usr/share/xml/udunits/udunits2.xmllibrary(ggmap)## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details.##
## Attaching package: 'ggmap'## The following object is masked from 'package:terra':
##
## insetlibrary(cartogram)##
## Attaching package: 'cartogram'## The following object is masked from 'package:terra':
##
## cartogramlibrary(patchwork)##
## Attaching package: 'patchwork'## The following object is masked from 'package:terra':
##
## arealibrary(tmap)
library(viridis)
library(tigris)## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.library(ggspatial)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggmap':
##
## wind## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutStatic Map 1: Location of attacks and population density
zombie.attacks.df <- st_as_sf(zombie.attacks.df)
?tm_polygons
tm_shape(nz.reg.shp.join) +
tm_polygons(col = "Region", palette = "viridis") +
tm_legend(outside = TRUE) +
tm_bubbles(size = "Population", col = "pop_density", border.col = "black")## Some legend labels were too wide. These labels have been resized to 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## The legend is too narrow to place all symbol sizes.
# required libraries
library(leaflet, quietly = T, warn.conflicts = F)
# start basemap
map <- leaflet() %>%
# add ocean basemap
addProviderTiles(providers$Esri.OceanBasemap) %>%
# add another layer with place names
addProviderTiles(providers$Hydda.RoadsAndLabels, group = 'Place names') %>%
# add graticules from a NOAA webserver
addWMSTiles(
"https://gis.ngdc.noaa.gov/arcgis/services/graticule/MapServer/WMSServer/",
layers = c("1-degree grid", "5-degree grid"),
options = WMSTileOptions(format = "image/png8", transparent = TRUE),
attribution = NULL,group = 'Graticules') %>%
# focus map in a certain area / zoom level
setView(lng = 174, lat = -41, zoom = 5)%>%
# add layers control
addLayersControl(overlayGroups = c('Place names',
'Graticules',
'Points',
'Lines',
'Polygons'),
options = layersControlOptions(collapsed = FALSE),
position = 'topright') %>%
# list groups to hide on startup
hideGroup(c('Place names'))
# show map
map### Add data
nz.shp.reproj <- st_transform(nz.reg.shp.join, crs = st_crs(zombie.shp))
library(htmltools)##
## Attaching package: 'htmltools'## The following object is masked from 'package:pander':
##
## plabs <- as.list(zombie.shp$Attacks)
labs2 <- as.list(nz.shp.reproj$pop_density)
## Color setup
colpal <- colorNumeric(palette = "magma", domain=nz.shp.reproj[['pop_density']], n=10)
colorData <- nz.shp.reproj[["pop_density"]]
map <- map %>%
addTiles() %>%
addPolygons(data = nz.shp.reproj, label = lapply(labs2, HTML), color = ~colpal(colorData))
map <- map %>%
addCircleMarkers(data = zombie.shp,
weight = 0.5,
label = lapply(labs, HTML),
col = 'black',
fillColor = 'darkslategrey',
radius = 4,
fillOpacity = 0.9,
stroke = T)
map